{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rosbag\n",
    "from cv_bridge import CvBridge\n",
    "import numpy as np\n",
    "import onnxruntime\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_images(bag):\n",
    "    #Read from bag\n",
    "    bag = rosbag.Bag(bag)\n",
    "    # Read Image and CompressedImage\n",
    "    bridge = CvBridge()\n",
    "    count = 0\n",
    "    imgs = []\n",
    "    for topic, msg, t in bag.read_messages():\n",
    "        if msg._type == \"sensor_msgs/Image\":\n",
    "            img = bridge.imgmsg_to_cv2(msg, desired_encoding='passthrough')\n",
    "            imgs.append(img)\n",
    "            count += 1\n",
    "        elif msg._type == \"sensor_msgs/CompressedImage\":\n",
    "            img = bridge.compressed_imgmsg_to_cv2(msg, desired_encoding='passthrough')\n",
    "            imgs.append(img)\n",
    "            count += 1\n",
    "    print(\"Read %d images from bag\" % count)\n",
    "    return imgs\n",
    "\n",
    "def compute_netvlad(sess, imgs):\n",
    "    # Compute NetVLAD with Onnx\n",
    "    netvlads = []\n",
    "    for img in imgs:\n",
    "        if len(img.shape) > 2 and img.shape[2] == 3:\n",
    "            img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n",
    "        else:\n",
    "            img_gray = img\n",
    "        img_gray = img_gray.astype(np.float32)[np.newaxis,:,:,np.newaxis]\n",
    "        output = sess.run(None, {sess.get_inputs()[0].name: img_gray})\n",
    "        netvlads.append(output[0].flatten())\n",
    "    netvlads = np.array(netvlads)\n",
    "    print(\"netvlads shape: \", netvlads.shape)\n",
    "    return netvlads\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read 47603 images from bag\n"
     ]
    }
   ],
   "source": [
    "bag_path = \"/home/xuhao/data/d2slam/drone3-RI/drone1.bag\"\n",
    "images = read_images(bag_path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "netvlads shape:  (47603, 4096)\n"
     ]
    }
   ],
   "source": [
    "netvlad_model = \"/home/xuhao/d2slam_ws/src/D2SLAM/models/mobilenetvlad_dyn_batch_size.onnx\"\n",
    "sess = onnxruntime.InferenceSession(netvlad_model, providers=[\"CUDAExecutionProvider\"])\n",
    "netvlads = compute_netvlad(sess, images)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "netvlads_pca shape:  (47603, 1024)\n"
     ]
    }
   ],
   "source": [
    "# PCA with sklearn of netvlads\n",
    "# Decompose netvlads to 1024 dimensions\n",
    "from sklearn.decomposition import PCA\n",
    "pca = PCA(n_components=1024)\n",
    "pca.fit(netvlads)\n",
    "netvlads_pca = pca.transform(netvlads)\n",
    "print(\"netvlads_pca shape: \", netvlads_pca.shape)\n",
    "# Normalize netvlads_pca\n",
    "netvlads_pca = netvlads_pca / np.linalg.norm(netvlads_pca, axis=1, keepdims=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dot_product err variance:  0.012703997\n",
      "dot_product = 0.647379 * dot_product_pca + 0.341641\n",
      "dot_product_pca = 1.461383 * dot_product + -0.499261\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAJRCAYAAAB1KtHzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABWNklEQVR4nO3dd7wcBdX/8c/ZvS29ElroBARRWkSaiIAKPiAIiKg0AQMK+khT/CGi8qBA1EdUFBGQItJRogTpiI/03kvogdBC+u13z++P2eyU5La9e3e2fN+vV16Z2ZndPZnce8+dmbPnmLsjIiJSLpm0AxARkfqixCMiImWlxCMiImWlxCMiImWlxCMiImWlxCMiImWVeuIxs4vM7F0ze6qX7WZmvzazOWb2hJltVe4YRUSkdFJPPMDFwO59bN8DmJb/MwP4fRliEhGRYZJ64nH3u4EP+thlb+BSD9wHjDez1csTnYiIlFpD2gEMwJrAG5H1ufnH5kV3MrMZBGdENDU2b736qqsV/YbZxiw9XT1FP79W6DgEdBwCOg6B+joOFlvLNYXnKm+8/PL77r5KMa9aDYlnQNz9fOB8gLE20ae9+bGiX+uAmXtw9Uk3lSq0qqXjENBxCOg4BOrpOGTGjImtPzdzk3DlqO++VvTrFh1R+bwJrBVZn5p/TEREqlA1JJ5ZwCH56rZtgUXuPq+/J4mISGVK/VKbmV0B7AxMNrO5wGlAI4C7nwfMBj4HzAFaga+lE6mIiJRC6onH3b/cz3YHjilTOCIidc0awrSwcK8Px7ZtdPR9heXXh/Ae1XCpTUREaogSj4iIlJUSj4iIlJUSj4iIlFXqxQUiIpKiTDa2mtsmLCgYe+WDw/OWw/KqIiIivVDiERGRQfmr3zCk5yvxiIjIgP3A72M0XUN6DSUeEREZkA18IZ9k7pBfR8UFIiJ1xhqbCsvZyRNj27rvfSJccY8tn8dtAPyI7YB7i35/nfGIiEi/fsVdAMxjJP+xNYf0Wko8IiLSp4/6u3yY+QAcxu5Dfj0lHhER6VXGc/yCuwE4gU+Ss6GnDSUeERHp1cX8E4CnmcQTVtSk6xWouEBEpMZFRx0AzDtmemF5tV/d0+vztvc3WZ1WAL7DziWLR2c8IiKygibv4cf5yrWj2Q3MSvbaSjwiIrKCG/gbAHcxlZdsfElfW4lHRERi9vCXaSD4DM8Ztm3JX1+JR0RECkZ6F8fzCACHlKB0emVUXCAiUoOi3Qky606NbVvtnN67DtxA0AD0r2zIPBs9LLHpjEdERAA40J8rLP/Othi291HiERERJnobR/AUAAew57C+lxKPiIhwFTcCcCGbscBahvW9lHhEROrcN/2xwvKV9qFhfz8VF4iI1IDs2LGx9dzG64TLjz4b3zky7mANX8oXmAPA3uw9fAFG6IxHRKSOXZLvxfZLtqbVGsvynko8IiJ16gd+HwBdGDfZemV7XyUeEZE6tIEvKIyx3pt9yvreSjwiIvXGnfO4HYDT2I4uy5b17VVcICJSpaLjDt46dLPYtjX+En4YtKe7O7YtOsb6niGOsS6GznhEROpIqcdYF0OJR0SkTgzHGOui4kjlXUVEpOyWl04/VcIx1sVQ4hERqQM7+puslh9jfVwJx1gXQ8UFIiJVwpqbY+tL99yisLzGbe/FtvXM/6Cw3OQ9nDZMY6yLoTMeEZEaN5xjrIuhxCMiUsOGe4x1MZR4RERqVDnGWBdDiUdEpEYtH2N9/TCOsS6GigtERCpYw5prFJbf/cw6sW2rPBAWEPQ8+2Js21c8HIXw+2EcY10MnfGIiNSYid7G13gaGP4x1sVQ4hERqTHlHGNdDCUeEZEaUu4x1sVQ4hERqRFpjLEuhooLREQqSKYlfmmsc8NVC8ur3Pd+bFuyoGB5L7afl3GMdTF0xiMiUgNO9aAlTicZbi7jGOtiKPGIiFS5DXwBO/EmAPtU8CW25ZR4RESqWcpjrIuhxCMiUsXSHmNdDBUXiIikLDruwNZfO7atYWF7YbnnuTmxbZUwxroYOuMREalC0THWx6c4xroY1ROpiIgUXOo3AfAkk3gyxTHWxVDiERGpMjv6XFbNj7E+PuUx1sVQ4hERqSJN3sNpXjljrIuh4gIRkRRYY1NhObPuWoXl7vEjYvtlHng6XHHnBv8rUDljrIuhMx4RkSpRiWOsi6HEIyJSBUZ5Z2GM9cFVVDq9MqknHjPb3cyeN7M5ZnbySravbWZ3mtmjZvaEmX0ujThFRNL01+7rAbiODXm7gsZYFyPVxGNmWeBcYA9gU+DLZrZpYrcfAFe7+5bAgcDvyhuliEi6vtwT3uc5r8LGWBcj7eKCbYA57v4ygJldCewNPBPZx4Gx+eVxwFtljVBEpASixQQAmfXCgoKOqeMKy03/eTq23/iuJXyNJ4HKHGNdDHP39N7cbH9gd3c/Mr9+MPBxdz82ss/qwC3ABGAUsJu7P7yS15oBzACYMG7i1medOrPouCZMHceCuYuKfn6t0HEI6DgEdBwCRR+HZNlzc5iIvCls7GlL22K7ff34wwF44HP78fgulXOnYcaJRzzs7tOLeW7aZzwD8WXgYnf/hZltB1xmZpu5ey66k7ufD5wPMNYm+tUn3VT0Gx4wcw+G8vxaoeMQ0HEI6DgEij0OfZ7xrDW+sBw94/lG54OF5VNuMripNo5/2sUFbwJrRdan5h+LOgK4GsDd7wVagMlliU5EJCVr5JawT8/zQGWPsS5G2onnQWCama1nZk0ExQOzEvu8DuwKYGabECSe98oapYhImf2p4wag8sdYFyPVS23u3m1mxwI3A1ngInd/2sx+Ajzk7rOAE4A/mtlxBIUGh3maN6ZERHqTuI9jDWHCyE5dPbatdf2JheUR9zxfWO5pb6+qMdbFSP0ej7vPBmYnHvthZPkZYIdyxyUikoZqG2NdjLQvtYmIyHKRMdY/ZPuqGGNdDCUeEZEK8b/L/gHAW4ziXlsj5WiGjxKPiEgF+Gj3PDbpCeqmvlblvdj6k/o9HhGRWmHZ+KWxaEFB60bxKaEjHn61sOyLFnI2wWd0gjHW1TdjZzB0xiMikrJL+SdQnWOsi6HEIyKSoh3aX67qMdbFUOIREUlJk3dz6uKbATiqSsdYF0OJR0QkJX997wIA7mQtXq7SMdbFUHGBiMgQWEP4YzQ7NV4CvWzTKYXlUc+8G9v2mfceIJsfY/1T+/gwRlh5dMYjIlJmo3raOS4/3eVg2yPlaMpPiUdEpMyufe1cAK5jWtWPsS6GEo+ISBkduOC+wvJ5mS3SCyRFSjwiImUyoXsphy74DwBftL1SjiY9Ki4QEelPtMzZLDZNNDslnEvZtmGiO8HcpYXl7ldf5y+5awC4wD7CQpqhTie86IxHRKQMjsk9Wli+yj6UYiTpU+IRERlma3R8wD7MAWBv2yfdYCqAEo+IyDD707NBFdvPbXrNjbEuhhKPiMgwOvWV4L5OrY6xLoaKC0REkhI906wpLCagsSFWULDso2sWlpvnt8eet8FT97Nj93NAfox1nRYTJOmMR0RkOLjzu+5bADgtu2PNjrEuhhKPiMgwOOeZoAHoW4zi3szUlKOpLEo8IiIltvmiV9h06VwADm/4r5SjqTxKPCIiJZTxHn7+3MUAnJDdhZzpx2ySigtERAAy4T0Ya4z/aMxOnlRYzo1oihUUZDtz4fMef4HL2v8KwJO2Ck/kJkKue7girlpKxSIiJbJjz+tMyY+xPqHp0ylHU7mUeERESqCpp4sfdv0bgKObPlc3Y6yLocQjIlICs+88DYA7M+vwcmZCytFUNiUeEZEh2nNuOMb6Z007phxN5VNxgYjUp2R3gkhBQXbC+Ni2pVtGPodjkGsOnzv5309ywpKgoOBgdsc7Okofa43RGY+IyBBcv+QvAFzX9OG6HGNdDCUeEZEiHfrcbYXl81u2STGS6qLEIyJShElti/j6s0Evti+NPjDlaKqLEo+ISBFu+OcZAFzYvDULMyNSjqa6qLhAROpHpKAg09wc3zRubGG5dfO1Yts6xoe/o2eWdXDSzZcV1q/sWh+6lpY60pqmMx4RkUEY++7b7N36JAD7Tjky5WiqkxKPiMggfOmM7wHwi7GfojXT1M/esjJKPCIiA3TWXRcD0EGWW0dukm4wVUyJR0RkADaeP5dPvfEUAPuvekTK0VQ3FReISO3KxMdNZ5oaw+VJE2PbWjdbo7C8ZGpjbNuqt77BZa/9CoBbDjuW9kvfCTe6lyjY+qEzHhGRfvzy7asAmNcwjtc22zLlaKqfEo+ISB+mvzWHTTrmAXDkmoelG0yNUOIREelFNtfD+bN/D8CJq31RY6xLREdRRKQXf78q6E7wZPOaPN0ytZ+9ZaBUXCAi1S864iByVpJpSXQnGBN2j+7ccNXYtkXrhwUFq93xHjssfo7Vli0C4PiuHeCNt/JbN1NBwRDpjEdEJKEp18Wpb/0NgG+se7jGWJeYEo+ISMJfX/glAHeO3ZRXWqakHE3tUeIREYnY7+l7C2Osz1rj8ylHU5uUeERE8sa0t3Lqv64F4ND1j045mtql4gIRqX69FBRkJoyP7da1bnjZ7IONW2LbVv3PAm5+/H8AuDazMfPeWAgsDDbmekoZbd3TGY+ICPCVd/5dWD4/q+4Ew0mJR0Tq3irLFnHo2/8C4EubfifdYOqAEo+I1L1b/vwTAC5YfRcWNo7uZ28ZKiUeEalrJ//7usLyNVO2TzGS+qHiAhGpDr10JwDIjAgLBWzkyMJy91qTY/t98KERheXJjy1lzfb5fOmZewDYp/kA/PlXAPDu7pKFLSvSGY+I1K2LH/81AD9ffx9arbGfvaVUlHhEpC6d9sKVAHRYAzevoiq2clLiEZG686H357LjgmcB+ML0k1OOpv6knnjMbHcze97M5pjZSr8CzOwAM3vGzJ42s7+UO0YRqSHu/OVv/wvADzf6Ml0ZXWIrt1SLC8wsC5wLfBqYCzxoZrPc/ZnIPtOA7wM7uPsCM1PHPpF6kOwI3de4g9GjCsu5NVYpLM//yKjYfhOfaeM3j50HwJs2hvvfyJJ548XgeR0dJQlb+pf2Gc82wBx3f9ndO4Ergb0T+3wdONfdFwC4+7tljlFEasQWC19ik6VzATiyea+Uo6lf5ikONDKz/YHd3f3I/PrBwMfd/djIPn8DXgB2ALLAj9z9nyt5rRnADIAJ4yZufdapM4uOa8LUcSyYu6jo59cKHYeAjkOgoo5DJvE7czay3hheOuseFT5uPT0cffhXAZh1wmm8s+Y68dfIDexnYUUdhxTNOPGIh919ejHPrYbP8TQA04CdganA3Wb2EXdfGN3J3c8HzgcYaxP96pNuKvoND5i5B0N5fq3QcQjoOATKfhwGc6ktMlk0eqnt/eljC8uzLw+6Ezw+dl1+fX8j2YfviL1Grr19QGHp62Ho0r7U9iawVmR9av6xqLnALHfvcvdXCM5+ppUpPhGpAbu8/ERhjPUJHzki5Wgk7TOeB4FpZrYeQcI5EPhKYp+/AV8G/mRmk4GNgJfLGaSIlEnkLMey2fim5si4gzHxfmq5yRMKy4s3HlNYHv9iB809nfz87ksAOHrEXmQfnxM8R8UEqUn1jMfdu4FjgZuBZ4Gr3f1pM/uJmS0f/XczMN/MngHuBE5y9/npRCwi1eamu38EwG2rbs4r2YmpxiKBtM94cPfZwOzEYz+MLDtwfP6PiMiA7fXm/YXln276JRoffD7FaGS5tO/xiIgMizGdrRz3wg0AfGXbE1OORqKUeESkJt16XXDh5Oq1duTtEbrEVklSv9QmInUmWSYd3dQQfgbHmuKtbKLjDnJTJsS2LdpkfGF55DudHPRKWCp90eKNaHokKCjoaRtYybQML53xiEhNmdSxmMNfuQ2AL435csrRyMoo8YhITbnmP2cC8IcNdmdRZkQ/e0salHhEpGac+ND1heWr1tkpxUikL0o8IlIT1lryHvvPCcZY77nTD/vZW9Kk4gIRSU20mADAIj3YMqPjIw18Qth3bekG42LbRs9t55r/OwuAX4z6JN1PzaMJyCWLCXI9JYhahkpnPCJS9U57NpgP2ZZp5NbmjVKORvqjxCMiVW3jBXP5xPxgduS+256ScjQyEEo8IlK1zHNcdOc5APxgk4PoyujuQTVQ4hGRqvWHu34LwNyWSdw36UMpRyMDpV8PRKT0ot0JLP77bXTcQWZES3zbyPBzNz5+TGxb67phQcGIt9vZYtHLfHjBGwDMGLM/DS8Eo7xyS5eFr9GjYoJKpDMeEak6Ge/h50//CYDvbHYkOdOPsmqi/y0RqTqXP/xLAB4bux5PjV0n5WhksJR4RKSqfHLeE6zSuRiAEz/8tZSjkWIo8YhI1Wju6eT0Ry4H4OubH9Nnp2upXCouEJGhSyaAyD2XTKQbAQCR4gIbMzq2yUeHow9a1xsf29a0qItb/+9UAO5onsYb83I08FbwvCVLIi+SC5fVqaAi6YxHRKrCnvMeKCyfPW63FCORoVLiEZGKN7qrle+8NAuAg6Yfn3I0MlRKPCJS8Wbf8mMArl5zR95u0RjraqfEIyIV7dAXbyssn7/e7ilGIqWi4gIRKU6koCDajQAgu/pqhWUfES8usNZwVIGPGRnb1rZ22J2gaWEnkzoWc8QLtwJw4JSv0fjSPAByra2x53lPLrKsgoJKpzMeEalYVz00E4Dz1/ksi7Ij+9lbqoUSj4hUpG+99PfC8tVTd0wxEik1JR4RqThTl73H3m8H5dOf/7hm7NQaJR4RqTiX//sXAJw1bV9aG1r62VuqjYoLRGTgeikoyEyYENvNR4bJwpa1xbdFCgo61hgb29a0qJPTnrsCgHZr4K72tWh8/f3gee2RooTOrvhrdkfW3Qf0T5H06IxHRCrGtKVvseOCZwE4YM2jUo5GhosSj4hUBPMcv3vyPABO3fgrdJkuyNQqJR4RqQjnPvh7ID/GeqLGWNcyJR4RSd1WH8xhk8VzAThii2+lHI0MN53LikjvEuMOrKmpsJwZH3YZ6F5/9dh+2WUdhWVvHBPb1j0+LDxoXNRBNtfDLx65EIATp3yRzFsLyQDelihKaO+IrOSIb1RBQTXRGY+IpOqyJ34FwGNj1uXpljXTDUbKQolHRFLziQ+eZnJXMMTtuxsfmnI0Ui5KPCKSiuaeTk596RoAZnz4GxpjXUeUeEQkFbPv+QkAt0/8CK+OXDXlaKScVFwgUu+SZxoW/j4aG3dgRna1KYXV1w9Yq7C82gPxQgCWhYvRYgKAbFsX//XuQ4X1n0/4LNkPlgLgnZ2Fx1foThAdfdDd3cs/RqqBznhEpKxGd7fx36/fCMCha3895WgkDUo8IlJW1z92NgDXrLod7zaO62dvqUVKPCJSNge/dkdh+Y9rfSbFSCRNSjwiUhaTOhZz2OtB4jlg8xNTjkbSpOICkXpn8d8/rTH8sRDtTkBDlnd3nVpYHfdKT2E52xq/2Z8bGXY4yC4LigSufiS4xHb+5J1Z3GpkaIW29tjz6Ih0PIgUGkBi9IFUNZ3xiMiw+1a+mADg+gkfSzESqQRKPCIyrNZsn89e7wfl0/tsfnLK0UglUOIRkWH1p2d+C8DMdfamNduccjRSCZR4RGTYFMZYZxq5ddIW6QYjFUPFBSL1IDneINKRwBriPwZsxIhwZfKEcLmhgZHvhgUFlouMIuiOjymwbpjWGo6x/uIG3yazLCgksPZod4JEAUFk9IF3JboTaPRBzdAZj4iUnLlz7nN/BOCHa+5HV0a/40pIiUdESu7Xz18AwBvNk3hg9IYpRyOVRolHREpqi8Uvs3HrWwB8fdNvphyNVCIlHhEpmaz3cPacywA4bqOvkTP9iJEV6cKrSK3KRAoIMonigkhBgbXES5xt1MjCci5RlJBrDNdHzgu7DmQ6g0KAy547B4DHRq7NszaFzLJ2rCNRQBBZjxYTQGLcQa4HqU36dURESuITi55lcncwV+fkqQemHI1UMiUeERmy5lwXp7xxPQBHbfh1jbGWPinxiMiQ3fBM0AD0jnGb8VrLlH72lnqnxCMiQ7Ln2w8Uls9ea+8UI5FqkXpxgZntDpwDZIEL3P3MXvbbD7gW+Ji7P7SyfUTqTvSSVnK8QbQ7QTbxO2ZjY7gt2qkA8BHNkeVwPzLGyDdbw9WO7mCM9Sv/AODQ9Y+GfLcB644UBiSLC9rDogTv6KO4QGpWqmc8ZpYFzgX2ADYFvmxmm65kvzHAfwP3lzdCEenLdU/9HIBrJmzDOxpjLQOU9qW2bYA57v6yu3cCVwIrO1c/HTgLaF/JNhFJwVffvruwfOGUT6UYiVQb8xQb75nZ/sDu7n5kfv1g4OPufmxkn62AU9x9PzO7CzhxZZfazGwGMANgwriJW5916syi45owdRwL5i4q+vm1QschULXHIXYZLrkt8jtnJvH7Z/SyXGR5wqRmFswPLo2NXPgBX/1+8G162Vm/o33U2PhrRH+u5BI/Y3K5lS8nn1ehqvbrocRmnHjEw+4+vZjnpn6Ppy9mlgF+CRzW377ufj5wPsBYm+hXn3RT0e97wMw9GMrza4WOQ6Cij0MJ7vFkRo+KbfLRkQ+Qjg+3ffFr07j64hcBuPX+0wA4f43duO6KeVjba/H3jt7jaW2LbcstC+8TeVt8WzXc46nor4cqkXbieRNYK7I+Nf/YcmOAzYC7LPgGWw2YZWafV4GB1KXk52MiySaaaIL1yLZEAYE1NxWWfVSiuCCyLXZG4pDp7OHYN2YXHrp+/Mewzu4Vz1QiIw18heKCyOiDHnUnqEdp3+N5EJhmZuuZWRNwIDBr+UZ3X+Tuk919XXdfF7gPUNIRSUl0jPUXNj4+5WikWqWaeNy9GzgWuBl4Frja3Z82s5+Y2efTjE1EVnTRs+cCMHPtvWnTGGspUtqX2nD32cDsxGM/7GXfncsRk4is6NPn/QKA1kwTt03aHOvoSjkiqVapJx4RqXzTlr3Fus8+DMABHzkx5Wik2inxiFS6gVauNca/na0p0p0gWjAA+JhIJVtkPwCy4ftZdy4YY/1sfoz1Ol+iK5cJig4iZdLWGT/7iVarrVi5Ftm3CsqnpfTSLi4QkQp3zgsXArBwyurcP3ZaytFILVDiEZFebbHkFTZumwfAtT88O+VopFYo8YjISmW9h7Ne+jMAx214GJ7scCBSJH0lichKXfb0rwF4dPS6PDN6rX72Fhk4FReIVII+CghiuyULCKLFBSNa4ttaIuvJAoLI63hD4v3c2WnhM0xaPsZ6g4PC14x0JIgWFES7EUC8W8EKbXBUUFD3dMYjIjHNuS5OeS0YYz1j46M0xlpKTolHRGJmPXkWALdN2IzXRmiMtZSeEo+IFOz5ftgGceba+6QXiNQ0JR4RAWB0dxvfmhu0+z9kk2P72VukeCouEElLL/dOLBN/3Boi36aJkubotuTog2hBgTfHiwtiBQU9wc3+5WOsr56yHe80TwhesysytsCJFQZ4ZM6OL10Wf/0OjT6Q3umMR0T46juRMdZr7JZiJFIPlHhE6tykrsUc8m6QeA7Y5LiUo5F6oMQjUuf+8lzwQdHzV9uNRQ2j+tlbZOiUeETq2LH5YgKA61bZNsVIpJ6ouECkXJLFBJEOBbEOBNnE74ONvY83iBUUJLoaeLQ7QeT1l5vaPp+95gczdvbe/OTCPpbLhTt1R4sLHFsWKSiIFhAkxyL0xJ8nEqUzHpF65M6Fz/8egLPX2UdjrKWslHhE6tCPXr0GCMdYi5STEo9InZnW+hbbLX4BgAM+fHzK0Ug9UuIRqSPmzm9fugiAU9f7El0Z3eaV8tNXnchwykSKBjJ9FBdECgFshREGYUGBjYyPPiBSNOAj4vdpPFrMkA2Wz3k+GGP9RvMk7p/woeA1u+KdBWKjDyLjDch5vFtBZBSCd8eLC1RQIH3RGY9IndhyycuFMdZf3/SbKUcj9UyJR6QOZL2HM1+6HIDjph1Gro9hcyLDTV99InUgNsZ6lMZYS7qUeERqXG9jrEXSouICkaGK3sRPXMKKFhTExhtArDDAmiIFBC2JD3NGxxskCg880snAG+PdCcw9GGP9Rn6M9Ye/gTcFMVhnpIAg2qkAINJ1wFvbw8dzuVi3gly0W4GKCWQQdMYjUsNmPf4zAG6dtDmvjlw15WhEAko8IjVqz/ciY6zX/0KKkYjEKfGI1KAx3W18a+5sAA7Z9NspRyMSN6h7PGb2X8CHgcKn2Nz9J6UOSkSG5ronzgbg6inb807z+HSDEUkYcOIxs/OAkcCngAuA/YEHhikukcrV13iDRHeCaEHBCkUD0Y4E0XEHyQKChkjRQENivEF2xe4EB715V+GhC9b9bP5Fer/5b+2dsXVf1hout0eKC9zxSFcDPFGUIDJAg7nUtr27HwIscPcfA9sBGw1PWCJSjEmdiznkrbsA+OKW3003GJFeDCbxLG/S1GpmawBdwOqlD0lEinXF478E4A9rfYZFjRpjLZVpMPd4/mFm44GZwCOAE1xyE5EK8O1X/15Yvm617VOMRKRvA0487n56fvE6M/sH0OLui4YnLBEZjKnt77Pne/kx1lt9P+VoRPo2mOKCY4DL3X2hu3eY2Ugz+6a7/24Y4xNJxwoFBJEOBNnEDf7obonuBDZiRLiSKBqIdiuIFg14NnEFPDIyIdmdwBsywRjr5/JjrNeLjLGOFBSsMPqgLVJQ0BEvLiDSnSA6+oBcLj7+QN0KpEiDucfzdXdfuHzF3RcAXy95RCIyKD+acyUAyzJN3DZ5i3SDERmAwSSerFn4a5+ZZYGmPvYXkWE2bdmbbL/weQAOUBWbVInBFBf8E7jKzP6QXz8q/5iIpMDcOfeZPwJwyrSvaIy1VI3BfKV+jyDZfCO/fiuqahNJza+fD779Xm+ZzIPj9ZE6qR6DqWrLmdmFwD1ADnje3Xv6eZpIdUqON4gWFCS3RW7+W7KAINqRoDG+jUgRQXTcgbckOhdECgo83xlhy0UvsdHyMdYfPQY3wxL3+mOjD5LFBZFxB75sWWxbbllbuC1ZTKCCAimBwVS1/RdwHvASYMB6ZnaUu980XMGJyIqyuR7Ofu5SAP570yM0xlqqzmAutf0C+JS7zwEwsw2AGwElHpEyuvyxoDvBI2PX55kxa6ccjcjgDeZXpSXLk07ey8CSEscjIn3Yaf5TTOoKxlh/70OHpByNSHEGc8bzkJnNBq4maJfzReBBM9sXwN2vH4b4RCSvuaeTU+dcA8DXP/LNFT/kKlIlBpN4WoB3gE/m198DRgB7ESQiJR6RYfT3R34KwC2TNcZaqttgqtq+1td2M/u+u/9s6CGJDKO+zhL6ukkfqWqLtboh0SanKVm5FqlIa0l83jpaydYQvnduROI1gL3eDkdfzdxg3/C9I0VmK1SuRdZtSWtsW3TOTq6tPb6tJ/I6qmKTYVDKcpgvlvC1RCRvTFcr337lHwB8dYvjUo5GZOhKmXh0wVlkGFz/0JkAXLnGjryrMdZSA0qZeHROLlJiB79xR2H5wnU+k2IkIqWjMx6RCjWpYzGHzL0LgP2nfy/dYERKaDCdC3Zw9//08dg1JY1MpBRWmKuTiSz2/rtSrGAgY7GCguTMnWiRgCXb4kTa6URn7gDkRq68ufvytjhXPvJzAH6/7u4sbBodbMzGY461xenoim9ri8zVaU0UF0QLCnoSna9y6oQlw2swZzy/6esxd//p0MMREYDvzLmhsHzdGjukGIlI6fV7xmNm2wHbA6uY2fGRTWOB3kcxikhRpra+x57vPATAXh//QcrRiJTeQC61NQGj8/uOiTy+GNh/OIISqVvuXPzorwE4c9p+tDU0q2xHak6/icfd/wX8y8wudvfXyhCTSN368YtXALAs28xtU7ZINxiRYTKYljkXmNkX3X0hgJlNAK50988OJQAz2x04h+Cy3QXufmZi+/HAkUA3QZuew5UApU8D7WGWjV8pjkx2j3cnsExszs4KBQSRmTvenCgYiMzcySU7F0TkGjJsvGQu2+XHWO+/zcnh+0W6B1hH/MZ/JlJQYEvbYtu8PVJc0BkvPPDOznC5uxuRchpMccHk5UkHwN0XAFOG8uZmlgXOBfYANgW+bGabJnZ7FJju7h8FrgXOHsp7ilQi8xznPhlMlf9/mxysMdZS0waTeHJmVhj+YWbrMPSrz9sAc9z9ZXfvBK4E9o7u4O53uvvyWtD7gKlDfE+RinPuE0HSeW3EKjwwUWOspbaZD7AJYP6S2PnAvwg+LPoJYIa731z0m5vtD+zu7kfm1w8GPu7ux/ay/2+Bt939f1aybQYwA2DCuIlbn3XqzGLDYsLUcSyYu6jo59eKmj8OK3zGJ7oc/k42YY3RLJgXGQ+9wvMi6318NohM4ve8/K5rPvskn/t10F/3j7/7M57YLzbSOvntGv3+7cnFt+VyK9+vv229qPmvhwHScQjMOPGIh919ejHPHUx36n+a2VbAtvmHvuPu7xfzpsUws4OA6YRjGZLxnU+QGBlrE/3qk4ofjHrAzD0YyvNrRdUehwF2oLbG+Jd/b/d49vvRTlx3xj2R55XoHk9Dhmyuh5seCpLOtzf7Os/8+ZUVeoBYboD3eFboQB25x9Oe6EDdFt4PGug9nqr9eigxHYehG0zngp3yi4vzf29qZrj73UN4/zeBtSLrU/OPJd97N+AU4JPu3pHcLhJLNtHkkiggiJ6RJLdFCwhihQdm8WSTSDze0hwuN/fxLdUQP5PxrHH5o8EY64fHbcAz4/JXshMnIJn2SGLojp/VxLoTRAoGAOiIJp74t40nuxWIlNFg7mCeFFluIbg/8zCwyxDe/0FgmpmtR5BwDgS+Et3BzLYE/kBwSe7dIbyXSEX55PtPMnH5GOsPH5pyNCLlM5hLbXtF181sLeBXQ3lzd+82s2OBmwnKqS9y96fN7CfAQ+4+C5hJ8AHWa/KXQl53988P5X1F0tbS08kPXgzaGx65+bEaYy11ZSg1m3OBTYYagLvPBmYnHvthZHm3ob6HSKWZ9fAZANyyyha8NnJIn0oQqTqDucfzG8KrzxlgC+CRYYhJpKbt9U5kjPWG+/axp0htGswZz0OR5W7giuSYBJFhNcBqtZhEeXN0pMEKhQeNkdEH0UKDjMUr15oSxQWR6jhPvKY3ZiLLWcZ0tfKt124E4CtbnYDn/00WKRqwZHVzpPTZEkUCRDsStMUr13KRde+Ody4YaAm1yHAYzD2eS4YzEJF68Nf7g+khV079BO+2jE83GJGUDGQswpP00aEg38pGRPpxyGu3F5YvWPezkNNZh9SngZzx7Jn/+5j835fl/z4INWwXGZDJHYs45I07Adjv499PORqRdA1kLMJrAGb2aXffMrLpe2b2CHDyyp8pIstd8VjwQdHfr7cHixpHpRyNSLoGU1xgZrbD8oICM9uewTUZFRmcQRQTWLSIIHKDf4XuBNFtI0fEXzNSeOAt8bEIPirc1xvjr5mLdCvwZHeCjHHci38rrF+3xvbhjf1os4VIn7UV2uIsifSJ645v8yVLwzg6Et0JogUFKiaQCjKYxHMEcJGZjSP4llkAHD4sUYnUiKmt77Hn2xpjLRI1mKq2h4HN84kHd1d7VpG+uHPJw+cA8LON9w/GWIvIwC+Vmdk4M/slcDtwu5n9YnkSEpEVnf7s5QAs1RhrkZjB3KO5CFgCHJD/sxj403AEJVLtNl4yl+0WBGOs99tWVWwiUYO5x7OBu+8XWf+xmT1W4nik3vUy3iBY7b3YINqRIFZc0BSfgWPReTmJ2TnRDgQkOhfkRkVGH2TjceSaw/dzs2CMdX6i6Pc/cgjd+TEKniiWyHRExh1Eph1kWuMdCKID3nzZstimaEGBdyXm6qigQCrUYM542sxsx+UrZrYD0NbH/iJ16fcP/Q6AV0dO4YFJG6ccjUjlGcwZz9HApZH7OgsADRERidj6gxfZaOlbABy5zbdTjkakMg0o8ZhZFjjY3Tc3s7EA7r64n6eJ1JVsroeZjwe3Pb+11VHkLIPpcpfICgaUeNy9Z/llNiUckZW76t6zAHhowoY8PW6dlKMRqVyDudT2qJnNAq4BCnc43f36kkcltW2A0zZXKCaIFBtYNnF7sjEsBogXGvS+nyeLC0aGBQS5pkh3gozRMzLSnSARf67B+NTbTzChK/i2OHHrwwv/RsuF+2YTN/+jxQXWHuky0NEZjytSUJBLjD7wzsi+OruSKjGYxNMCzAd2iTzmgBKP1LWWnk5Oe/IKAA7b7jsaYy3Sj8F0LvjacAYiUq3+ecdpwd+rb8Wro1dNORqRyjeYzgXrm9nfzew9M3vXzG4ws/WGMziRSrf33HsLy2du9sUUIxGpHoP5HM9fgKuB1YE1CO71XDkcQYlUgzFdrfz3i38H4Es7fjflaESqx2Du8Yx098si6382s5NKHZDUiIHe5+ijO0GsSABiHQnIJJ4XLRSIdi4Y0RLbz0dEOhA0xV8/WlAQG3Vghkfi6mkK3vuGO/8HgMvX+STvNk/AIt0HCmF2hQ9m2uLFBdGCAlsWfhbb25MFBJHCg574WAQVFEg1GkziucnMTiY4y3HgS8BsM5sI4O4fDEN8IhXpay/eWli+YMPdU4xEpPoMJvEckP/7qMTjBxIkovVLEpFIhVulfRGHvXQ7AJ/fRTN2RAZrMFVtfRYS5Edj39rXPiK14Nq7fgbAbz60J4uaRpPp0uUukcEo5ejqs0r4WiIV6bjn/lpYvnbdHfvYU0R6M5hLbf3Rp+bq2WA+NBnrQJCNb4t2K0gUEETXLdF1IFaI0Fd3gsZIAcGIxti2npbwNXIN0TigpznL2kvfZa+3HgTgM587nVxDPh4Lz3iy7fGb/5mucD3T2hHbZktbw7giBQXeR3cCTxYXiFShUp7x6HqD1C53Lvv3LwA4fcsDNcZaZAhKmXhEatbPHr4EgCUNLdyy1lYpRyNS3QbTuWCFX/ESj71aioBEKs0qr77E9u89C8A+u56acjQi1W8wZzz39vWYu+879HBEKot5ji/8NCiZ/u70r9GdKeVtUZH61O93kZmtBqwJjDCzLQmLCMYCI4cxNqli0aIBz8Vv/1lj+GVnyaKEaEFBovAgWkBgjfHCAB81IlyJvH6uJV5ckBsZPi/XFP+9q7slEnO+uODCf/0GgFfGrMo9UzcJQuyO/3uyHeEN/4ZlXbFtmaVhoYC1JooGOsJiA28Pl3Pt8SIEPNISQZ0KpAYM5Ne3zwKHAVOBX0YeXwL8v2GISaQiTH/3BTZaFIyxPmSX41KORqR29Jt43P0S4BIz28/drytDTCKpy+Z6+N/7LgDgb9//Cf6s6nBESmUw3023m9kvzeyh/J9fmNm4YYtMJEXX33IGAA+sMo13Ntgo5WhEastgEs+FBJfXDsj/WQz8aTiCEknTLm89xsTOpQCcsO2RKUcjUnsGU6KzgbvvF1n/sZk9VuJ4pJL11Z3Aev8dZoXuBFHJAoKmSDFAJv5+1hIZcdCY+NKNdCjIRTsQtMSLEKIFBT2J4gIMWro7+dFjwRjrgz95fNCUwIn15WhojXcPaFgaFhREiwkgXlDgra2xbdGCgmihATl1J5DaNpgznjYzKzSnMrMdgLY+9hepOrfcHHxOZ/bUrXltjMZYiwyHwZzxHA1cGrmvswA4tPQhiaRj31fvKSyfucUBfewpIkMxkM/xHB9ZvRQYlV9eBuwGPDEMcYmU1ZjOVr7z9A0A7L/LySlHI1LbBnLGMyb/98bAx4AbCK54HwQ8MExxiZTVP+74CQB/3mBn3h0xIeVoRGrbQD7H82MAM7sb2Mrdl+TXfwTcOKzRSTqiRQQDLCiwRCFAfPRB4lZipKAgVkwAsYKCWDEBQFNk3EGiuMCz4fNyzeG2nuZ48UKuMdyva3Sw7Xd3/q7w2Hkf/Vz+BSNPMmhYGt7wzya7E7SGYwusvTO2zaPFBYmOBN4Zvo7GHUg9GUxxwapA9LuqM/+YSNVad/E7bDH/FQD2++z3U45GpD4MprjgUuABM1s+gnEf4OJSByRSTn+55ecA3DZ1c94ZqUtsIuUw4MTj7meY2U3AJ/IPfc3dHx2esESG36x//KSw/MNtDyLTowacIuUwqB7v7v4I8MgwxSJSNh9Z8CqT25cA8NnP/zjlaETqi4aLyIoFBL203l+hA0G0W0GyuCC6W7KAIFoYYH10J2hKjD4YEc4d9MZ4LD2RcQceKWZYPt5gue6RWXDnt/88D4A/T9uZZdkRZHoci9zfz3RGRhHknMZId4Ls0niRgC0LP0ftydEHkW4F0WICSBQUaNyB1BG13JW685+/frewfN5mn0sxEpH6pMQjdWWnN58sLO/y+TNSjESkfinxSN3IeI6f3X8pAP+7+d50Zhv7eYaIDAclHqkbd9wSDsy9doMd+9hTRIaTigvqVV8dCTKRzgINfZwVRLsMNCS+lDKR32mSIwwi722jRsa3RQoYvCVelJAbEcYS7U4A8Q4Fnok+Hqzs/eq9hcd23O8s3CATbzIQG3fQEOlOkOlxskvCggJrS3QgiI43aE8UF/TVnUDjD6RO6YxHal5TTxcnPhF87vmUbQ8m18fsIBEZfvoOlJp3+42nFJbvnPrRFCMREVDikRp35LP/LCx/Yq+zUoxERJZT4pGaNbqrjUNfvAOAb21/VN/3tUSkbFRcUC+SP3QHMtLALD7SINP77ykrdDWIdB2w5ubEzuH7eXOigGBUpHNBNh5XT0v45Zprir9f98hwvac5eN4/rg16sS1qGsnDa2wY/BO6ww4Bma54t4BsW3e4rS3SZcA9VlDgS1ujT4sVFKww+iBaUKBiAhGgAs54zGx3M3vezOaY2QqjH82s2cyuym+/38zWTSFMqTKnPnBFYXkP9WITqSipJh4zywLnAnsAmwJfNrNNE7sdASxw9w2B/wV0oV76tErrQvZ4Pehle9Cnj+9nbxEpt7TPeLYB5rj7y+7eCVwJ7J3YZ2/gkvzytcCuZrpYL727YXbQCuepiWvz8rjVU45GRJLMU+yKa2b7A7u7+5H59YOBj7v7sZF9nsrvMze//lJ+n/cTrzUDmAEwYdzErc86dWbRcU2YOo4FcxcV/fyq00sen7DmWBa8tTi648BfI7reR+fqFe4bRdeTTbNjI7kT2/Lv8fmzTmP1Oc8DcN6FV64YZvTLPfGlbz1hR2rLhRsnTG5hwTuR+zqR/YLXibxQro9tVa7uvi96oeMQmHHiEQ+7+/RinlszxQXufj5wPsBYm+hXn3RT0a91wMw9GMrzK0a0A0Hih79HfrAmCwOWFxTsf8YuXHvKHeGGxngXg1i3gqaBb/PIem50vPAg1xJu6xqT6JoQ+SHeMyKRsBzWW/w2R+WTzr67/z/eufGNWBcDgKZFYQFBtiOeJBoi3Qkyi8NEs+93Nue6M8LOB97aFnued3audDkZc7Wrme+LIdJxGLq0L7W9CawVWZ+af2yl+5hZAzAOmF+W6KSqXHbHLwG4ZeoWGmMtUsHSTjwPAtPMbD0zawIOBGYl9pkFHJpf3h+4w9O8PigV6W83nV5Y/vE2X00xEhHpT6qX2ty928yOBW4GssBF7v60mf0EeMjdZwEXApeZ2RzgA4LkJFLw0fdfYXJHMMZ6j8/9KN1gRKRfqd/jcffZwOzEYz+MLLcDXyx3XFIl3Pndv38PBGOslzSN7OcJIpK21BOPlNAgqswz0Rv+ycqy5cUGlomNKYh1MQCsJVIYkOhAEC0g8JZ4kUBuZLhv96j4tp6m8D1yTfF/T64hE9kv2Hb/X04sPPbbrfcMwu8Mr8Q2Lo13C8h2hQUF2Wh3AsCWRUYatEWWcx5bTxYQxNZ1FVikX2nf4xEp2iffCMdYf3L/n6UYiYgMhhKPVKVMLsfZ/w4+Vzxz633o6GtgnYhUFCUeqUr3XvndwvK1G2uMtUg1UeKRqrPvS/cUlrc78OwUIxGRYqi4oNr1UVAQ61aQHFsQXU8UF0SLCKIFBDYyXjHmzZECgpHxDgTeGL5+17iW2LZo14FoMQFArjGMuWtEsg0PNHV3ceKjwRjrk3c6mJ7G4PnRVjjRdjfZzkR3gkVhdwJrjY8wYPHSMI5od4Jcjly0uKAnMd5ABQUig6IzHqkq91wWTs64Y53NU4xERIqlxCNV45sPhx/32uag4pvAiki6lHikKozuaOPwJ24HYMYe39QYa5EqpsQjVeGuK34AwMLmkTyy2gYpRyMiQ6HigmrTVzFBHwUEK4w+aIz81zckvgwsU3ivaEGBj4wXCeTGhOu5pvjrR0cadI6Lv35P5CM3nWOT3QnC9WxHcNP+9Dv+Unhs14NOj21bLhO539/QFhYUNCxNdCfoDMci2JJlsW0eLSBojxQe5HLxggJPzNwRkUHRGY9UtClLF7Lniw8DcMAXTuxnbxGpBko8UtFuvjw4w3liyjq8NFFjrEVqgRKPVKyLbvhNYfnQL3w7xUhEpJSUeKQibTD/bbZ8+1UAdv/qD9INRkRKSsUFlShZQGCR3w8SN7atj+aYFi0aSI4+iGyzxsRrjMgXDWQzsYKC7smjY7v1tES6E4yJFxe0jwvXOybE/z2dY8PlbKJ5QEO+YcBfLw9a4dy0wZa8O2oC5mCR+/vZrnhxQbaXgoLMwngBgbWHIwx8WWtsW6/dCdwhl+hWICJF0xmPVJw7zy/MAeSUXQ9KMRIRGQ5KPFJRtnrzZSa1BT3TPnno6SlHIyLDQYlHKoc7F1/7WwAumL4LS5o1xlqkFinxSMV45qwTCsu/3mHPFCMRkeGk4oJK0VfvsUhBwQrdCaIv0Rj/74zuayNHxHeOFBckOxKw/HUyGbonjSo83DkuXoTQOTZ8/bZJ8d9hWlcPb/53jUvcmM+F/9bm+cHzdnv2icJj0791Jp6JFxMANLZFxh20xYssGheHBQUNi8KRBrY0XkAQ7UjgnZ3xbdGCAhUTiAwbnfFI6jK5HOdefTEAZ3xqXzoamtINSESGlRKPpO7Z08NWOFduoTHWIrVOiUdS9ZUH/1NY3uTUn6cYiYiUixKPpKa5q4vTZl8HwLe+eCi55IdcRaQmqbggTb0VFNggfgBnwtewpvjNfxsVFgbQEC9K8FFhsUGuOf68XEvwZeFZo3N8eL+lfWL8NZatEca5dJN4C4KN1n6nsPzGgvGxbe1vBXE9f8L3C4/duvHmQTFBpCFBsrggE+lW0NAa35htC4sLbElYUJAsIMi1RrZ1dRPfqIICkXLQr5iSihP/flNh+UOn/CLFSESk3JR4pOzGtrZy7C13AHDQwcdojLVInVHikbJ74nunATB/9CgeWkdjrEXqjRKPlNXZ51xXWN76Zz9KLxARSY2KC6pArANBNvG7QmSkgbX00oEA8BHNsU09o8P17jHxD2x2jwjez7MZ2iaFr9G2SvyS2LJ1wpvxq622MLbtlfcmFpa7lgWvv/oHC/n83U8CsMdxJ5JtzWDxBgSx9Yb2+OiDpsVhMUC2NV40YEvDbgXeFllOjD6IFRQkRkyISHnojEfK5r4TzwTgkbXX4YXVNMZapF4p8UhZXPfT3xeWv3iMxliL1DMlHhl20958h+lzXgNg25+fnHI0IpI2JR4Zdred+r8A/G3bLZg3cXy6wYhI6lRckKZeOhRYJn4TPzbuIDEWIdatoCH+3xktKMiNihcXdI0P16PjDQA6RwVx5RqgddUwlo6J8Zv9ngnX33l5cnxbNtj28Ik/Ljz2na8chC0FIv88T3wFRgsKGlvjN/+z7WExQybSqQDi4w9yrZHigu5Ed4JoQYHH/z0iUh4645Fh87EXX2HykmUAfPTMn6QcjYhUCiUeGR7uXPvzoKDg3E/vwuKRGmMtIgElHhkWrx39vcLyzL0+l2IkIlJplHik5HZ/NBxjvdFvz0gxEhGpRCouKKdMojAgM8CxCJEmmsnOBdYUdh3wkfHOBbmWcFuyO0G0oKBjTPw1e/Iv4xnojlwhyzXHb8ZbdxiXNwU37TO5HOf98TIATj1obzpGZaEz0SGgI3zvhmXxTdnIdIVM4nnZZWG3AlvWFtvm0YKCSHcCzyUKCFRQIJI6nfFISb1ydDhj59Jdt0sxEhGpVEo8UjKH3HlPYXm9C3+aYiQiUsmUeKQkmju7OP2KGwA4+qiDNMZaRHqlnw5SEi8c+4PC8k1bfyTFSESk0qm4YLgNdLpmdPRBU7wQwKKv0RjfRkvYgcBHxrsTdE4Kiw26xsb/q9vHhb9z9LTEY2yflH+9RqNnRHgzPllc4CODTgInXfvPwmMfuvY0RlgbbYvD97ae+OtnIhMNsu3xf07zorA7QcOyeNeB2OiDpfGqBO8MX9R7wtcg14OIVBad8ciQjF3WyrE33gXAQT85TGOsRaRfSjwyJE8eG7TCmT92JA99eN10gxGRqqDEI0X71R+uLCzv8KfvphiJiFQTJR4pyuoLFvKF+x4DYLf/OS7dYESkqqi4oNSS9zh6GX2Q3BYtILDka0RGH1iiO4FHuxOMjhcedI8KCxY6R8fj6BrT+7iDsHOB0zU+vDnfMjm8uX/vMUErnCc2XYP2bRroWZj4d0a6GmRa49uaF4TbWhbGuxM0LokUFyyOVx7ERh+0d8S2xcYfeKJTgohUFJ3xyKBd8f8uKCx/81dfTTESEalGSjwyKNNef4etXngDgH3/fFTK0YhINVLikUH5x/G/A2DWJz7Ku1PGphyNiFQjJR4ZsHsOP7uwfNJ/75diJCJSzVRcUAoD/NBkcgxCbMRBpHMByXEJkW0rjD4YGS0uaIxt6xoVvn7bKvHX7JgQFhR4fFoD3ZO6gscbnMlrLQRgy6deZ9LioFvA564/llXHLOGD1hGF5/T0JMY1dIbrjcvi7924JHzvpiXJ4oKwA0FmcWtsm3dEuhN0x7saxLoVaPSBSEXTGY/0z50LvvdnAC798sdZOqalnyeIiPQutcRjZhPN7FYzezH/94SV7LOFmd1rZk+b2RNm9qU0Yq13D+/5s8LyBYd/IsVIRKQWpHnGczJwu7tPA27Prye1Aoe4+4eB3YFfmdn48oUou/znucLybjd+J71ARKRmpJl49gYuyS9fAuyT3MHdX3D3F/PLbwHvAquUK8B6Zz05Zv70egDO/MZn6WzSLUERGTrzlG7EmtlCdx+fXzZgwfL1XvbfhiBBfdh9xY+mm9kMYAbAhHETtz7r1JlFxzZh6jgWzF1U9PMTgQ1sW2/LANGhatn47woeWc81xp/nkTyRa4hvyzWG/+/WkLjB3xDcqD/q8wcXHrvoH5dgxL9WOnrCN+jqilcoWLRzQVdiLEKXr3QZwLoiRQI9iZEGPZE4k1+3w/h1XNKvhyqm4xDQcQjMOPGIh919ejHPHdZfYc3sNmC1lWw6Jbri7m5mvf7kMLPVgcuAQ1eWdPKvcT5wPsBYm+hXn3RT0XEfMHMPin5+JvEDOFqhlk1si65H5/E0J+bxjAirx3zsqNi2rokjC8utq8fn8bRPCN+7bUpi5s4qkdY0U9pi29YYs5h9//ZIYf2Au79OLnsPLdmu2H4vLZpcWH7r/YmxbdkFYYXdyHnx9x79ZvhfOOrteOubxreXhCsLF8e2+ZKl4XJk/g4Mb1XbkL4eaoiOQ0DHYeiGNfG4+269bTOzd8xsdXefl08s7/ay31jgRuAUd79vmEKViKaOLr5z7h0AzDzj0+SyKn4UkdJJ8yfKLODQ/PKhwA3JHcysCfgrcKm7X1vG2OrabXueU1i+f+f1U4xERGpRmonnTODTZvYisFt+HTObbmbLu1AeAOwEHGZmj+X/bJFKtHXixL/cWlj+ww2XphiJiNSq1MqU3H0+sOtKHn8IODK//Gfgz2UOrX+DGO/suchN/Gzv+8VGISTuBdEQrnePjX94s2tM+F/YMTYeV9foyD2eNeOf9F9rvfcKywetfT8ATYu6OPyGewH49xUbsGrzEraf8FJhvwcWrhd7jfau8L2tIx5z8wfhezctit9zaVoa3o/JLo7fq7HWcBRCrjV+7yl6X8eThQfqViBSNXTxXgoO/3iQdNonNTD/Y6P62VtEpDhKPALAric8W1j+5/2bphiJiNQ6JR5h1Lx2pt0YXHq76u9bpxyNiNQ6JR7h4E89AMDbW45lwTRdYhOR4aUeKMVI3sjuo9jAkoUCUdFtzeGHP60l/kHQ3NjwQ6K55vjrdY6NFB6MSH5INIxzw2nzYttOXnc2AB/b77XCY1NmT2EGb/FSV/BBzZcz3cxpXbWw/ZXF8Q+JfvDm+MJyy7vxuBqXhcsjFiQ6IywOCx0ySxMFBK2RUQhd8Q+sRgs1VEwgUr10xlPHRj3fwfhHgyqy3MPrphuMiNQNJZ46tv0erwLg+42BNXTyKyLlocRTpz659ZzCsv921T72FBEpLSWeOvSRx+fStCD4AOadj26YcjQiUm90faWMrK+OB5nexyJ4Q/j7Qfeo+E38aLeC1tXjN9xza4c37ldpyXd2dudXx10DwKNHTWXOiNVYveul2PMuX/hxADbpHsW/3wx7tS2eH694a/wgjKVpYfyfM+L9sKCgYVm8y0DjB5ECgsVLY9u8Lexc4N3xbgvkEt0KRKQq6Yynzly5wwWF5YeOWze9QESkbinx1JFt7nylsHzR49unGImI1DNdaqsT1pPj+FNuA+DCE3Yg1zzwRqciIqWkM546ccUnLiws37qferGJSHp0xlNisVHXALFJ3b3neWsI/ys8Mfq6Z0Q4RrorUVzQPil8v56R8Q4By7ccfHM4uHXnfx5PblmG29mk8NjN8z8ce95bS8cBsGZXC4s/CAsKmuY1xvZrmR++d/PCxOiDxWEsLe+0xrbZ4rCtQbSYABKjD5LFBSJSE3TGU+OaO7v48cX/AOCU0/bWGGsRSZ1+CtW4Zw/9cWH53ztOSzESEZGAEk8N++4VNxeWd7rlhBQjEREJKfHUqHHLWjl61r8BOOC0Iwc1rltEZDipuKDEYq37WUmxQXRbb/dbGhKjFPpIGg2R+/bZtnC/x/7fDwF4b8xo7p86jez78RED7y8NiwaWLWqJbcssCooIulZroPmNsNCheUH8vVs+CP+tI9+NFwI0LewIX29xvLjAl4bFBbmOjvi2HnUnEKl1OuOpQYf+392F5em/PDXFSEREVqTEU2MmLFvKqf+4AYBtz/p+ytGIiKxIiafGPHj6aQD8ZpdPM2/i+HSDERFZCSWeGnLCP2cXls/5zO4pRiIi0jsVF5Sax7sHYI0r3w8gUYhQ0J24we7hfpnu+HMalwbrqy1dwDfuuh2AnY48nVFvGt0LwvfOZeNxdI8IX2dEZ7x4IZu/35+ZZIycF+7XvDj+3s0LwjibPoh3IGh4f0kY/qLFsW251rDYwLsS3Qm8l2MiIjVDZzw1YvaV/wPAmTt9gUUjRvWzt4hIepR4asBpd19ZWL5i80+kGImISP+UeKrcegveZu8XHgRg26N/lnI0IiL90z2eaubOddfPBODkTx1EW2NzygGJiPRPiWe4RYsNMn0c7j4+sZ/tCLc1Lwg7EJz52J8AWNg0iv+M/iij34q/Rueo8ITWE80QEkHG1hrag/Xsh5xR74TxNy5JdCdYEBYUZBYluhMsiXQniHQqgERBQU6dCkTqjS61ValNF73ONh+8CMAXPn1KytGIiAycEk8VMs/x24fPA+CEjx9Bd19nUiIiFUaJpwr98YHfADBn9Oo8MGXjlKMRERkcJZ4q87H5L7D+sncAOOpjx6QcjYjI4OkaTalZ77k82fLfGhoi28Kb+NYa7wKQyY9PaMj1cNazFwPwrU2OJLu0m+bItIPuUfH/zmxbpLAhMVkh1xQ+kOmIFxdkO4PnZbqdEe+EsTQsbIvHvyxc92QBwbJIsUGycEIFBSJ1TWc8VeSK588B4P5x03h+9NSUoxERKY4ST5XYZeFTjMkFZx+nTvtKytGIiBRPiacKtPR08t03ZwFwxGbHaIy1iFQ1JZ4q8Lfnfg7A7Alb8MaIVVKORkRkaFRcUGrJsQj03jIgWmxg0Rvu3eEn+/de+nhh+TcTdiXT2hl7DesM9820x98r1xT+92a6E3FFxw/0xIsLMvnXtK4cjW8vCje0xYsevL0jfK/WROeCaHeCFY6JiNQznfFUsLE9bRy9+N8AHLz2jJSjEREpDSWeCnbVOxcC8Jfx2/J+w5iUoxERKQ0lngp16OJ7C8uXTdwhxUhEREpLiacCrdK9hAOXPgzAAasdkXI0IiKlpeKCUojeqE+UOnsu3GZ9jSaIFBpc+v4lAJw7egcW9zRg3ZEihLaO+POy4e8OlnjvTB9l19YZaXng8eICOvIFDN09+IKwuMA74u8dLSBIdmVQdwIR6Y3OeCrMccvuLiz/fcRmKUYiIjI8lHgqyDo9C/hsZzBjZ+/Jh6ccjYjI8FDiqRTu/GHx9QD8bMyudFhjygGJiAwPJZ4KccbSmwFYaC38q2XDlKMRERk+Ki4ohWJ7p+VvyG/S8x5bd78JwEHjDoSe+Cf9vTPsVmDd8d8VLBupWEje4I8UDXjyNbs6e91GV77woKebXGTcwQoFBNGOBMkCBRGRXuiMJ2Xmzq86g7Od74/ene4+S99ERKqfEk/Kzuv4BwBzshN5tHHNlKMRERl+Sjwp+ljPm6zrwedkvjVmn3SDEREpEyWelDR4D//TeScA327eHdeMHRGpEyouKIU+OhfERG7OX9V1HQD3Z9bg+cxkiHYSyMbv81jkJr5b4neF5A3/aFh9bYu8X6/75Rzv7qPDgYhIEXTGk4Jdc68yhuAH+g+bPpVyNCIi5aXEU2Yt3sX3eu4D4MjmvTTGWkTqTmqJx8wmmtmtZvZi/u8Jfew71szmmtlvyxnjcJjVHVxim23r80ZmXMrRiIiUX5pnPCcDt7v7NOD2/HpvTgfu7mN7Vdi35/nC8q8atkkxEhGR9KRZXLA3sHN++RLgLuB7yZ3MbGtgVeCfwPQyxVY6+cKAcd7B0f4oAF/O7BWMS4h2DIgtJ272D/By3AodCGLbiuw6oIICESmxNM94VnX3efnltwmSS4yZZYBfACeWM7DhcK3PAuDPtinv28iUoxERSY/5MP5Ga2a3AautZNMpwCXuPj6y7wJ3j93nMbNjgZHufraZHQZMd/dje3mvGcAMgAnjJm591qkzi457wtRxLJi7qP8dB2j67OvY8o7ZAPzx5xeGGwZaWDDQ+oO+/iuL+H8u9XGoVjoOAR2HgI5DYMaJRzzs7kVdhRrWxNPnG5s9D+zs7vPMbHXgLnffOLHP5cAngBwwGmgCfufufd0PYqxN9I/brkXHdsDMPbj6pJuKe3IimazirfzFbwRg38w+LLHmcNdsL33ZMvHXSE4W7U2pL7UN6TjUEB2HgI5DQMchcJtfW3TiSfMezyzgUODM/N83JHdw968uX46c8fSZdCrN8qTzW9sylnREROpVmonnTOBqMzsCeA04AMDMpgNHu/uRKcZWvMjZw0n+YGH5Bltxxk6vHQOSD0fPgJKdC6JdDfroVKAiARGpFKklHnefD6xwPczdHwJWSDrufjFw8bAHViLr+CI+w2sA7GVfSDkaEZHKoc4Fw8GdC7gVgDPs47SbWuKJiCynxDMMfsb/AbCAZu6ytVOORkSksijxlNim/j7TeQeAL/NfKUcjIlJ5dA2ohDLunMNdAHzPdqInOcbaey93jkkUEMSLBlRAICLVTWc8JfSH/H2dOYznEVuhEYOIiKDEUzLb+DzWZTEA37TdUo5GRKRyKfGUQIPnOIP/AHCM7aox1iIifVDiKYGr+TsA97I6L9jElKMREalsKi4Yol39tXCMNduvuEMxN/y9jwICEZEqpzOeIRjhXZxM0BbncD6jMdYiIgOgxDMEs/J9TW9kPd6wsSlHIyJSHZR4irSvv1BY/pVtnWIkIiLVRYmnCOO8g2/wBAAHqjuBiMigqLigCNfmq9guYxPm24j4RnUPEBHpk854Bulwf7KwfKl9OMVIRESqkxLPIEzxZXyZ5wHYl8+nHI2ISHVS4hmEywnmrP+WLVhiTSlHIyJSnZR4Bqi/MdYiIjIwSjwDEBtjzT7pBiMiUuWUePoTGWP9P2iMtYjIUCnx9ONM/g0EY6z/ZWulHI2ISPVT4unDpv4+W/MuoDHWIiKlosTTi9gYaz5Bj+lQiYiUgn6a9mL5GOsXNcZaRKSklHhWYq1nnyiMsT6GXVOORkSktijxJDR6D7tfeA4Ax7CLxliLiJSYEk/CNRpjLSIyrJR4Inbz1xhFN9DLGGsRERkyJZ68Ed7F9/JjrK/+7hkaYy0iMkyUePKWj7H+B+uzaMpqKUcjIlK7lHiIj7E+x7ZKMRIRkdpX94lHY6xFRMqr7hNPn2OsRUSk5Gqy1fJGW6/PrQ9d0/+OJ58MZwWLB/szHJx/+K677uLW3ACeX+N0HAI6DgEdh4COQ8CGUIBVv2c8r70GZ+Wzzvz56cYiIlJH6jfxrLtu8Pe558JEfVBURKRc6jPxHHZYuPzNb6YWhohIPaq/xPPUU3DJJcHy0qXpxiIiUofqK/G4w0c+EixfeSWMGpVuPCIidai+Es9nPhP8PWUKfOlL6cYiIlKn6ifx3HMP3HZbsDx3brqxiIjUsfpIPD09sMMOwfKtt0JjY7rxiIjUsfpIPMvv62y1Fey2W7qxiIjUudpPPDfeCM8+Gyw/+GC6sYiISI0nno4O2HPPYPmBByBT2/9cEZFqUNs/iadMCf7eay/42MfSjUVERIBaTjyXXQaLFwfLN9yQbiwiIlJQm4knl4NDDgmWn3tOY6xFRCpIbSaeRx8N/j7qKNh443RjERGRmNpMPMudd17aEYiISEJtJp7Ro+GNN9KOQkREVqI2E8/GG8PUqWlHISIiK1GbiUdERCqWEo+IiJSVEo+IiJSVEo+IiJSVEo+IiJRVaonHzCaa2a1m9mL+7wm97Le2md1iZs+a2TNmtm6ZQxURkRJK84znZOB2d58G3J5fX5lLgZnuvgmwDfBumeITEZFhkGbi2Ru4JL98CbBPcgcz2xRocPdbAdx9qbu3li1CEREpuTQTz6ruPi+//Daw6kr22QhYaGbXm9mjZjbTzLLlC1FERErN3H34XtzsNmC1lWw6BbjE3cdH9l3g7rH7PGa2P3AhsCXwOnAVMNvdL1zJe80AZgCsuuqqW1955ZVFx7106VJGjx5d9PNrhY5DQMchoOMQ0HEIfOpTn3rY3acX89yGUgcT5e679bbNzN4xs9XdfZ6Zrc7K793MBR5z95fzz/kbsC1BMkq+1/nA+QDTp0/3nXfeuei477rrLoby/Fqh4xDQcQjoOAR0HIYuzUtts4BD88uHAiub1vYgMN7MVsmv7wI8U4bYRERkmKSZeM4EPm1mLwK75dcxs+lmdgGAu/cAJwK3m9mTgAF/TCleEREpgWG91NYXd58P7LqSxx8Cjoys3wp8tIyhiYjIMFLnAhERKSslHhERKSslHhERKSslHhERKSslHhERKSslHhERKSslHhERKSslHhERKSslHhERKSslHhERKathHYuQFjN7D3htCC8xGXi/ROFUMx2HgI5DQMchoOMQ2NjdxxTzxNR6tQ0nd1+l/716Z2YPFTtnopboOAR0HAI6DgEdh4CZPVTsc3WpTUREykqJR0REykqJZ+XOTzuACqHjENBxCOg4BHQcAkUfh5osLhARkcqlMx4RESkrJR4RESmruk48Zra7mT1vZnPM7OSVbG82s6vy2+83s3VTCHPYDeA4HG9mz5jZE2Z2u5mtk0acw62/4xDZbz8zczOryZLagRwHMzsg/zXxtJn9pdwxlsMAvi/WNrM7zezR/PfG59KIcziZ2UVm9q6ZPdXLdjOzX+eP0RNmttWAXtjd6/IPkAVeAtYHmoDHgU0T+3wTOC+/fCBwVdpxp3QcPgWMzC9/o16PQ36/McDdwH3A9LTjTunrYRrwKDAhvz4l7bhTOg7nA9/IL28KvJp23MNwHHYCtgKe6mX754CbAAO2Be4fyOvW8xnPNsAcd3/Z3TuBK4G9E/vsDVySX74W2NXMrIwxlkO/x8Hd73T31vzqfcDUMsdYDgP5egA4HTgLaC9ncGU0kOPwdeBcd18A4O7vljnGchjIcXBgbH55HPBWGeMrC3e/G/igj132Bi71wH3AeDNbvb/XrefEsybwRmR9bv6xle7j7t3AImBSWaIrn4Ech6gjCH7DqTX9Hof8ZYS13P3GcgZWZgP5etgI2MjM/mNm95nZ7mWLrnwGchx+BBxkZnOB2cC3yhNaRRnszw+gRlvmyPAws4OA6cAn046l3MwsA/wSOCzlUCpBA8Hltp0Jzn7vNrOPuPvCNINKwZeBi939F2a2HXCZmW3m7rm0A6t09XzG8yawVmR9av6xle5jZg0Ep9PzyxJd+QzkOGBmuwGnAJ93944yxVZO/R2HMcBmwF1m9irB9exZNVhgMJCvh7nALHfvcvdXgBcIElEtGchxOAK4GsDd7wVaCBqI1pMB/fxIqufE8yAwzczWM7MmguKBWYl9ZgGH5pf3B+7w/B21GtLvcTCzLYE/ECSdWryeD/0cB3df5O6T3X1dd1+X4F7X59296EaJFWog3xd/IzjbwcwmE1x6e7mMMZbDQI7D68CuAGa2CUHiea+sUaZvFnBIvrptW2CRu8/r70l1e6nN3bvN7FjgZoIKlovc/Wkz+wnwkLvPAi4kOH2eQ3CD7cD0Ih4eAzwOM4HRwDX52orX3f3zqQU9DAZ4HGreAI/DzcBnzOwZoAc4yd1r6krAAI/DCcAfzew4gkKDw2rtF1Mzu4Lgl4zJ+XtZpwGNAO5+HsG9rc8Bc4BW4GsDet0aO04iIlLh6vlSm4iIpECJR0REykqJR0REykqJR0REykqJR0REykqJR0REykqJR6QfZvYjMzuxj+2HmdkaZYxn3d7a1A/guTub2faljklkMJR4RIbuMGDIiSfflmm47Qwo8UiqlHhEVsLMTjGzF8zs/4CN849tke/G/ISZ/dXMJpjZ/gSNUy83s8fMbEQvr/eqmZ1tZk+a2QNmtmH+8YvN7Dwzux84e2Xvkd9vazN73MweB46JvO5hZvbbyPo/zGzn/PLuZvZI/nm3WzDI8GjguHysnyj9kRPpnxKPSIKZbU3QHmkLgnYgH8tvuhT4nrt/FHgSOM3drwUeAr7q7lu4e1sfL73I3T8C/Bb4VeTxqcD27n78yt4jv8+fgG+5++YD/DesAvwR2C//nC+6+6vAecD/5mP990BeS6TUlHhEVvQJ4K/u3uruiwkaIY4Cxrv7v/L7XEIwnXEwroj8vV3k8WvcvcfMxq3sPcxsfP7xu/OPXzaA99oWuDvfPRp372uYl0hZKfGIlI/3srxsCK/ZTfz7uGUIryVSFko8Iiu6G9jHzEaY2RhgL4LksCByX+RgYPmZyRKCeT39+VLk73uTG9190creIz9gbaGZ7Zh//KuRp70KbGFmGTNbi2BkMwRjG3Yys/UAzGziIGMVGTZ1OxZBpDfu/oiZXQU8DrxLMJsFgtlM55nZSIL5M8tbwF+cf7wN2K6P+zwTzOwJoINgeuXK9PYeXwMuMjMHbons/x/gFeAZ4Fngkfy/4T0zmwFcn5+e+i7waeDvwLVmtjfBPSPd55Gy01gEkTLITy2d7u7vpx2LSNp0qU1ERMpKl9pESsjM/gqsl3j4e/lx2SKCLrWJiEiZ6VKbiIiUlRKPiIiUlRKPiIiUlRKPiIiU1f8H+fWN0SDMRPEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dot_product dot_product_pca\n",
      "0.5 0.23143065880125185\n",
      "0.6 0.3775689572494522\n",
      "0.7 0.5237072556976525\n",
      "0.8 0.669845554145853\n",
      "0.9 0.8159838525940535\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Verify the PCA result, we random choose two netvlad vector 1000 times\n",
    "# and compute the dot product before and after PCA by compute the variance\n",
    "import random\n",
    "dot_product = []\n",
    "dot_product_pca = []\n",
    "for i in range(1000000):\n",
    "    idx1 = random.randint(0, len(netvlads)-1)\n",
    "    idx2 = random.randint(0, len(netvlads)-1)\n",
    "    dot_product.append(np.dot(netvlads[idx1,:], netvlads[idx2,:]))\n",
    "    dot_product_pca.append(np.dot(netvlads_pca[idx1,:], netvlads_pca[idx2,:]))\n",
    "dot_product = np.array(dot_product)\n",
    "dot_product_pca = np.array(dot_product_pca)\n",
    "print(\"dot_product err variance: \", np.var(dot_product-dot_product_pca))\n",
    "plt.figure(figsize=(10,10))\n",
    "heatmap, xedges, yedges = np.histogram2d(dot_product, dot_product_pca, bins=100)\n",
    "extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]\n",
    "plt.imshow(heatmap.T, extent=extent, origin='lower')\n",
    "plt.xlabel('dot_product')\n",
    "plt.ylabel('dot_product_pca')\n",
    "plt.grid(which='both')\n",
    "# Set to equal scale\n",
    "# Fit line of dot_product_pca to dot_product\n",
    "m, b = np.polyfit(dot_product_pca, dot_product, 1)\n",
    "print(\"dot_product = %f * dot_product_pca + %f\" % (m, b))\n",
    "m, b = np.polyfit(dot_product, dot_product_pca, 1)\n",
    "print(\"dot_product_pca = %f * dot_product + %f\" % (m, b))\n",
    "plt.plot(dot_product, m*dot_product + b, 'r')\n",
    "plt.show()\n",
    "#Print a table of original and PCA dot product\n",
    "print(\"dot_product dot_product_pca\")\n",
    "for i in range(5, 10):\n",
    "    thres = i / 10\n",
    "    print(thres, m*thres+b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pca_result shape:  (1025, 4096)\n"
     ]
    }
   ],
   "source": [
    "pca_result = pca.components_.copy()\n",
    "#Add the mean vector\n",
    "pca_result = np.vstack((pca.mean_, pca_result))\n",
    "print(\"pca_result shape: \", pca_result.shape)\n",
    "# Save the pca result as csv\n",
    "np.savetxt(\"../models/netvlad_pca.csv\", pca_result, delimiter=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pca_result_sp shape:  (65, 256)\n"
     ]
    }
   ],
   "source": [
    "# Load the pca result\n",
    "pca_comp_sp = np.loadtxt(\"../models/components_.csv\", delimiter=\",\")\n",
    "pca_mean_sp = np.loadtxt(\"../models/mean_.csv\", delimiter=\",\")\n",
    "# Save as stack \n",
    "pca_result_sp = np.vstack((pca_mean_sp, pca_comp_sp))\n",
    "print(\"pca_result_sp shape: \", pca_result_sp.shape)\n",
    "np.savetxt(\"../models/superpoint_pca.csv\", pca_result_sp, delimiter=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1.46 * 0.8 - "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.5 ('base')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "08ce52785f0fedc81003ce387e097a83d6cc9494681cd746006386992005bb71"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
